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SECTION  1.  INTRODUCTION 

Considerable  effort  has  been  invested  in  developing  methods  and  models  for 
predicting  acoustic  wave  propagation  in  an  ocean  environment  in  which  variabil¬ 
ity  in  only  two  dimensions,  range  and  depth,  is  considered.  Since  direction, 
speed,  and  intensity  of  an  acoustic  wave  in  seawater  are  affected  by  three- 
dimensional  environmental  conditions  such  as  seamounts,  fronts,  and  eddies, 
variability  must  be  accounted  for  in  the  third  dimension  (azimuth).  A  three- 
dimensional  parabolic  equation  (PE)  that  accounts  for  azimuthal  variability  in 

an  ocean  environment  was  developed  by  Tappert  during  the  last  decade.  Several 

2  3 

years  later,  Baer  and  Perkins  *  introduced  methods  that  use  the  fast  Fourier 
transform  to  implement  the  three-dimensional  PE.  Later  still,  Siegmann, 
Kriegsmann,  and  Lee^5^  (SKL)  developed  a  three-dimensional,  wide  angle  wave 
equation  of  which  the  three-dimensional  parabolic  approximation  introduced  by 
Tappert  is  a  special  case.  Mathematical  models  for  solving  the  SKL  three- 
dimensional,  wide  angle  wave  equation  were  then  developed  by  Schultz,  Lee,  and 
Jackson^  and  by  Lee  and  Siegmann? 

Because  existing  three-dimensional  models  are  limited  in  capabilities 
and/or  are  slow  due  to  the  nature  of  the  numerical  methods  employed,  the  search 
continues  for  an  accurate  and  efficient  three-dimensional  model  with  a  variety 
of  useful  capabilities. 

g 

A  mathematical  model  recently  developed  by  Lee,  Saad,  and  Schultz  (LSS) 
offers  an  accurate  and  efficient  approach  to  solving  a  new  version  of  the 
three-dimensional,  wide  angle  wave  equation.  Their  approach  is  based  on  the 
successful  application  of  an  Implicit  Finite  Difference  (IFD)^’  scheme  for 

solving  two-dimensional  acoustic  wave  propagation  problems.  Lee,  Saad,  and 
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Schultz  extended  the  utility  of  the  IFD  technique  to  include  three-dimensional 
effects . 

In  this  report,  we  develop  a  computer  model  (F0R3D)  that  implements  the 
LSS  method  for  solving  the  LSS  version  of  the  three-dimensional,  wide  angle 
wave  equation.  The  model  is  capable  of  predicting  acoustic  propagation  loss  in 
range-,  depth-,  and  azimuthal-dependent  ocean  environments.  An  important  fea¬ 
ture  of  the  model  is  that  it  is  also  capable  of  accurately  treating  wide  angle 
wave  propagation  in  the  vertical  plane.  Computational  speed  is  favorable  since 
the  LSS  method  requires  only  solving  two  tri-diagonal  systems  of  equations  for 
each  step  marched  forward  in  range. 

The  acronym  F0R3D  selected  for  the  model  represents  the  numerical  methods 
employed;  i.e,  finite  (F)  differences  methods,  ordinary  (0)  differential  equa¬ 
tions,  and  rational  (R)  function  approximations  to  solve  the  LSS  three-dimen¬ 
sional  (3D)  wide  angle  wave  equation. 

The  derivations  of  the  SKL  and  LSS  three-dimensional,  wide  angle  wave 
equations  are  presented  in  sections  2  and  3,  respectively.  Section  3  also 
includes  the  mathematical  formulation  of  the  LSS  method  for  solving  the  LSS 
three-dimensional,  wide  angle  wave  equation.  Section  4  presents  the  finite 
differences  methods  used  to  employ  the  LSS  method while  Section  5  summarizes 
the  algorithm  and  computations  programmed  in  the  computer  model.  The  structure 
of  the  computer  model,  including  a  brief  description  of  each  subroutine,  is 
given  in  Section  6.  Section  7  is  a  miniature  users’  guide  and  is  devoted  to 
input  and  output  formats.  A  test  problem  for  demonstrating  accuracy  and  the 
capabilities  of  the  model  is  included  in  Section  8. 
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The  model  is  written  in  Fortran  for  a  VAX  11/780  computer  and  has  been 
designed  such  that  additional  capabilities  may  easily  be  incorporated.  Prelim¬ 
inary  test  runs  have  been  made  on  other  computers,  such  as  the  VAX  8600/FPS  164 
and  the  CRAY  XMP  12. 

User  contributions  that  will  enhance  the  model  are  invited.  The  develop¬ 
ment  of  the  model  continues. 


3/4 

Reverse  Blank 


TR  7943 


SECTION  2.  SKL  THREE-DIMENSIONAL,  WIDE  ANGLE  WAVE  EQUATION 


The  formulation  of  the  SKL  three-dimensional,  wide  angle  wave  equation  is 
summarized  below.  A  complete  development  of  the  equation,  including  discussions 
on  approximations  and  conditions,  can  be  found  in  references  4  through  7. 

2.1  DERIVATION 

If  we  assume  a  geometry  which  is  cylindrically  symmetrical  as  shown  in 
figure  1, 


Figure  1.  Geometry  of  Cylindrical  Wedge 


the  three-dimensional  Helmholtz  equation  can  be  written  as 


2  2 

+  V  P  =  o. 


(2.1) 


where 


p  =  p(r,0,z)  =  acoustic  pressure 


n  =  n(r,0,z)  =  c^/c(r,0,z)  =  index  of  refraction 
Cq  =  reference  sound  speed 
c(r,0,z)  =  sound  speed  in  the  ocean 


=  w/c  =  reference  wave  number 
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0)  =  27Tf 

and 

f  =  source  frequency. 

Using  Tappertrs  parabolic  decomposition  technique,  we  let 

p(r,0,z)  =  u(r , 0 , z )  v(r ) ,  (2.2) 

in  which  the  factor  v(r)  represents  the  rapid  variation  of  the  pressure  that  is 
modulated  by  u(r,6,z).  Substituting  equation  (2.2)  into  equation  (2.1),  we 
find  that 


[v  +  —v  ]u  + 
rr  r  r 


[u  +  (—  +  )u  +  u  +  ^ru 
rr  r  v  r  r  zz  2  00 

r 


k^n^(r,0 ,z)u]v  =  0.  (2.3) 


Using  as  a  separation  constant,  and  setting  the  expression  in  the  first  set 

2 

of  square  brackets  in  equation  (2.3)  equal  to  -k^v  and  the  expression  in  the 

2 

second  set  of  square  brackets  equal  to  k^u,  we  obtain 


v  H - v 

rr  r  r 


+  v  -  o, 


(2.4) 


and 


u  +  {—  *t - 

rr  r  v 

Equation  (2.4)  is  the 
in  terms  of  out-going 
first  kind 


v  )  u  +  u  +  —  u  +  k  (n  (r,0,z)  -  l)u  =  0.  (2.5) 

r  r  zz  2  00  0 

r 

zero-order  Bessel  equation  for  which  the  solution  for  v 
waves  is  given  by  the  zeroth-order  Hankel  function  of  the 


v(r  ) 


■  Ho1>(kor) 


(2.6) 
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Introducing  the  far-field  approximation 

kQr  »  1,  (2.7) 

enables  us  to  replace  the  Hankel  function  by  its  asymptotic  expansion  in  equa¬ 
tion  (2.6),  giving 


v(r) 


=  Hq1} (kQr) 


i(kQr  -  tt/4) 
e 


(2.8) 


Using  equation  (2.8)  to  simplify  the  coefficient  (1/r  +  (2/v)  v^)  in  equation 
(2.5),  we  obtain  a  partial  differential  equation,  which  is 


u  + 
rr 


2iknu 
0  r 


+  u 

22 


ee 


+  k^(n2(r ,0 ,2) 


l)u  =  0. 


(2.9) 


Expressing  equation  (2.9)  in  operator  form,  we  find 


(fc  +  ik0  '  ik0  \A  +  (n2  -  1) 

\  k 


ai. 

2  3z2  (kQr)2  362/ 


(h  +  ik0  +  ^Vl  +  (n2  -  1)  + 


2  2 

a  +_i  a 


9  9  9 

0  3z  (kQr ) Z  30 


u  =  0.  (2.10) 


The  first  operator  in  equation  (2.10)  represents  the  out-going  wave,  and  the 
second  operator  represents  the  in-coming  wave.  Considering  only  the  out-going 
wave,  we  obtain 


(fc  +  iko  '  iko  yf-  +  (n2  '  1}  +7 


3  2 


(kor)‘ 


30‘ 


u  =  0. 


(2.11) 
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Defining  the  operators, 


x 


(r  ,9  ,z) 


(2.12) 


and 


1  3 


2  2 5 
(kQr)  30 


(2.13) 


the  three-dimensional  out-going  wave  equation  (2.11)  may  be  rewritten  as 


( §7  +  iko  ’  iko  /  1  +  x  +  y)U  =  °* 


(2.14) 


Using  a  rational  function  approximation  for  the  square  root  operator,  we  have 


/T+lT+y  =  y 


!  +  Plx  +  p2y 


+  +  q2ys 


(2.15) 


where  p^,  q^,  and  are  constants  which  influence  the  size  of  the  propaga- 

12 

tion  angle  in  the  vertical  plane.  Claerbout's  coefficients,  p^=P2=3/4  and 
q  =q  =1/4,  ensure  an  approximation  of  the  square  root  to  second  order  in  x  and 
y.  Substituting  equation  (2.15)  into  equation  (2.14),  we  have  the  SKL  wide 
angle,  three-dimensional  parabolic  equation, 


('ik0 


+  ik 


1  +  PjX  +  p2y 

01+  qxx  +  q2y 


(2.16) 


derived  by  Siegmann,  Kriegsmann,  and  Lee. 
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2.2.  SOLUTION 

Solutions  for  the  SKL  three-dimensional,  wide  angle  wave  equation  can  be 
found  in  references  4  through  7,  and  are  not  discussed  here  since  this  report 
is  devoted  to  solving  the  LSS  equation. 


9/10 
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SECTION  3.  LSS  THREE-DIMENSIONAL,  WIDE  ANGLE  WAVE  EQUATION 


In  the  previous  section  it  was  shown  how  Siegmann,  Kriegsmann,  and  Lee 
derived  their  equation  by  substituting  a  rational  function  approximation  for 
the  square  root  operator  in  the  one-way,  out-going  wave  equation  (2.14).  By 
using  a  higher  order  approximation  for  the  square  root  operator,  Lee,  Saad,  and 
Schultz  derived  another  version  of  the  three-dimensional,  wide  angle  wave  equa¬ 
tion  that  may  then  be  solved  by  successive  tri-diagonal  systems  in  the  result¬ 
ing  finite  difference  equations.  Their  derivation  and  method  of  solution  are 
summarized  below.  For  complete  details,  the  reader  is  referred  to  reference  8. 

3.1  DERIVATION 

Lee,  Saad,  and  Schultz  begin  their  development  with  the  three-dimensional 
out-going  wave  equation  (2.14), 


(I—  +  ikn  “  ikn  /  1  +  x  +  y)u  «  0 
3r  0  0 


(3.1) 


as  previously  derived  by  Siegmann,  Kriegsmann,  and  Lee.  It  is  at  this  point  in 
the  development  of  the  three-dimensional,  wide  angle  wave  equation  that  the 
approach  introduced  by  Lee,  Saad,  and  Schultz  varies  from  that  of  Siegmann, 
Kriegsmann,  and  Lee. 

Based  on  the  assumption  that  azimuthal  variation  is  small  but  not  negligi¬ 
ble,  Lee,  Saad,  and  Schultz  introduced  the  higher  order  approximation, 


/  1  +  x  +  y  =  (1  +  -|x  -  -|x2)  +  \y 


(3.2) 
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where  the  terms  inside  the  parentheses  on  the  right-hand  side  of  the  equation 
give  wide  angle  capability  in  the  vertical  (r,z)  plane.  Substituting  equation 
(3.2)  in  (3.1),  we  obtain  the  LSS  three-dimensional,  wide  angle  wave  equation, 

h  u  =  (_iko  +  iko  [1  +  ix  '  I*2  +  iy])u-  (3-3) 

Equation  (3.3)  is  the  equation  that  is  solved  by  the  computer  model  presented 
in  this  report. 

3 . 2  SOLUTION 

Locally,  the  solution  of  equation  (3.3)  can  be  expressed  analytically  by 
an  ordinary-differential  equation  solution, 

-ik  Ar  ik. ( 1  +  Lc  -  -^x2  +  ly)  Ar 

u(r  +  Ar , 6 , z )  =  e  e  u(r,6,z)  (3.4) 

Assuming  that  n(r,0,z)  varies  slowly  with  respect  to  0,  the  operators  x  and  y 
are  nearly  commutative,  allowing  equation  (3.4)  to  be  rewritten  in  the  form, 

-ikQAr  ikQ(l  +  -|x  -  -|x2)Ar  ikQ  -|yAr 
u(r  +  A r,0,z)  =  e  e  6  u(r,0,z)  (3.5) 

Setting  6  =  ikQAr  and  rewritinS  equation  3.5,  we  have  the  iterative  scheme, 

-6  6(1  +  4x  -  -jU2)  -y 

u(r  +  Ar)  =  e  e  e  e  u(r)  (3.6) 

developed  by  Lee,  Saad ,  and  Schultz,  in  which  the  exponentials  are  approximated 
by  rational  function  approximations. 

Substituting  the  approximations, 
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6(1  +  —x 


F2) 


6  1  +  +  f)x 

-  e  [ - j - f-] 

1  +  (4  '  4)x 


(3.7) 


and 


v 


1  +  f y 

i  -  *y 


(3.8) 


into  equation  (3.6)  results  in  the  final  expression. 


1  +  (y  +  |)x  1  +  fy 

u(r  +  Ar  )  =  [ - : - j — ]  [ - j—  ]  u(r) 


1  +  (}  -  f)x 


i  -  fy 


(3.9) 


The  advantage  of  using  the  approximations  (3.7)  and  (3.8)  is  that,  since  x 
and  y  are  real  and  small,  then  both, 


[1  +  <i  -  £)  x] 


(3.10) 


and 


11  -~y] 


(3.11) 


are  non-zero,  thus  assuring  that  both. 


6  1  +  (i  +  |)x 

5  l+(i-|)x 

4  4 


(3.12) 


and 


1  -v> 


(3.13) 


are  unitary.  This  property  assures  the  stability  of  equation  (3.9) 
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In  the  next  section,  finite  difference  methods  are  used  to  develop  a 
numerical  solution  for  equation  (3.9). 
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SECTION  4.  SOLUTION  TO  THE  LSS  EQUATION  BY  FINITE  DIFFERENCES 

In  the  previous  section,  the  derivation  of  the  LSS  three-dimensional,  wide 
angle  wave  equation  (3.3),  and  a  method  (equation  (3.9))  for  solving  the  LSS 
equation  were  presented.  In  this  section,  finite  difference  methods  are  em¬ 
ployed  to  obtain  a  numerical  solution  for  equation  (3.9).' 

4.1  AN  ALTERNATE  FORMULA 

Lee,  Saad,  and  Schultz  began  by  rewriting  equation  (3.9)  in  the  form, 


j+1  ri  ,  ,1  ,-1  ri  ,  A  ,  fi\  i  ri  $  ,-1  ri  .  i  j 

uJ  »  [1  +  (.J  -  ^)x]  [1  +  (4  +  ^)x]  [1  -  -^y]  [1  +  £ylu 


(4.1) 


where  represents  the  acoustic  field  u(r,0,z)  at  range  r,  u^  represents  the 
acoustic  field  u(r4-Ar,0,z)  at  range  r+Ar,  and,  as  previously  defined, 


6  =  ikQAr, 


(4.2) 


o  -j  3  2 

x  =  (n  (r  ,0  ,  z)  -  1)  H - -  - -  and 

ko  9z 


(4.3) 


y  = 


1  d< 


(kQr)2  902 


(4.4) 


By  exploiting  the  near  commutativity  of  x  and  y,  Lee,  Saad,  and  Schultz 
derived  an  alternate  formula  to  equation  (4.1), 


[1  +  (|-  -  f)x]  [1  -  |-y]uj+1  =  [1  +  (7-  +  |)x]  [1  +  |y]uj 


4  4y 


(4.5) 


that,  for  reasons  discussed  in  reference  8,  is  unconditionally  stable. 
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4.2  TWO-STEP  SOLUTION 

Given  the  acoustic  field  u^  at  the  present  range  r,  the  object  is  to  solve 
for  the  field  u**+\  one  step  forward,  at  the  advanced  range  r+Ar.  This  can  be 
accomplished  by  solving  equation  (4.5)  in  two  steps. 

Writing  RHS  (right-hand  side)  for 

[1  +  (^  +  |;)x]  [1  4-  |*y]u^ ,  (4.6) 

and 

vj+1  for  [1  -  ^y]u^+\  (4.7) 

the  two-step  marching  process  is  as  follows: 

Step  1 

4+1 

Compute  RHS  and  solve  for  v  in  the  following  system  of  equations: 

[1  +  -  f)x]vj+1  =  RHS.  (4.8) 


Step  2 


Now  that  v^+^  is  known,  solve  for  the  unknown  u^+^  in  the  following 


system  of  equations: 


[l  -  ^y]u^+1  =  (4.9) 

The  distinct  advantage  of  this  two-step  procedure  is  that  only  two  tri-diagonal 
systems  need  to  be  solved  for  each  step  marched  forward  in  range.  Consequent- 
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ly,  less  memory  and  fewer  computations  are  required  than  were  for  previous 
methods  . 

4.3  DISCRETIZATION  OF  OPERATORS 

Before  the  two-step  procedure  previously  explained  can  be  performed,  the 
operators  x  and  y  are  discretized  by  central  differences.  In  equation  (4.8)  of 
step  1,  the  x  operator  is  replaced  by  its  discrete  analogue  as  follows: 

Writing  equation  (4.8), 

[1  +  (i  -  f)x]vj+1  =  RHS ,  (4.10) 


2  1  O 

and  substituting  (n  (r,0,z)-l)  +— — -  for  the  operator  x,  we  have 

ko  32 


(1  +  (■£•  -  |o  [n2(r,e  ,z) 


1  +  \  ^-r]  )vj+1  =  RHS 

ko  32 


(4.11) 


32  v^+1 

Replacing  - t—  by  its  discrete  analogue, 

3z 


32  J+1  j+l+  J+l  , 

*  2  fK  v2  nrfl,£  m ,  £  m-l,£ 

9  z  (Az) 


(4.12) 


and  rewriting  equation  (4.11),  we  have 


(1  +  (7  -  7)  (n2(r,0  ,z  )-l))  v^+J  + 
4  4  t  in  m,£ 


+  (|  -  i)  -5-^ - j  (vj^j  -  2vi+l  +  vi+l  )  =  RHS 

4  4  ,2/a  x2  nrfl.fc  m,£  m-l,£ 

Kq(Az; 


(4.13) 


17 


TR  7943 


where  Az  is  the  depth  increment,  and  j,  m,  and  £,  are,  respectively,  in¬ 
dices  in  range,  depth,  and  azimuth.  Equation  (4.13)  may  then  be  expressed 
in  the  matrix  form, 


A^+1v^+J  =  RHS , 
m,  2  ’ 


(4.14) 


where  the  tri-diagonal  matrix  contains  the  diagonal  elements: 


uoDer  diagonal: 


k^(Az)2 


(4.15) 


main  diagonal:  1  +  (t-  “  tO  (n2  -  1)  -  2(v-  -  t)  x  ~ — j  (4.16) 

ko(4s) 


lower  diagonal: 


k2(Az)2 


(4.17) 


The  solution  vector  contains  L*M  elements  where  m  =  1,2,...,M  for  each  2  - 

m ,  Z 

1,2,...,L.  M  is  the  number  of  acoustic  sensors  in  the  vertical  column,  and  L 
is  the  number  of  vertical  columns  as  varied  from  port  to  starboard  in  azimuth. 

In  equation  (4.9)  of  step  2,  the  y  operator  is  replaced  by  its  discrete 
analogue  as  follows: 

Writing  equation  (4.9), 

[1  -  fy]uj+1  =  Vj+1  (4.18) 

1  3 2 

and  substituting  —  -  — for  the  operator  y,  we  have 
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,,  6,1  3  NN  j+1  j+1 

(1  -  t  < - J  — ))u  =  v  ' 

4  (k  r )l  902 


(4.19) 


92  uj+1 

Replacing  - - —  by  its  discrete  analogue, 

36 


a2  uJ+1  -  1  (uj+!„  -  2uj+i  +  uj+J  ,), 


30 


(A0) 


2  m,£+l  m,£  m,#,-!' 


(4.20) 


and  rewriting  equation  (4.19),  we  have 


j+1  _  6_ 
Jm,  SL  4 


(uj+1  -  2uj+1  +  uj+1  )  =  vj+1 

2  , .  „  x2  m,£+l  m,£  m,£-l  m,£ 


(kQr)  (A0 ) 


(4.21) 


where  A0  is  the  azimuthal  increment,  and  j,  m,  and  i  are,  respectively, 
indices  in  range,  depth,  and  azimuth.  Equation  (4.21)  may  then  be  ex¬ 
pressed  as  a  system  of  equations  in  the  matrix  form, 


i+1  j+1  j+1 

BJ  uJ  =  vJ  n 
m,  £  m,  £ 


(4.22) 


where  the  tri-diagonal  matrix  contains  the  diagonal  elements: 


5  1  1 

upper  diagonal:  -  ^  j  2  - j 

kQr  (AG) 


(4.23) 


main  diagonal:  1  + 


6  1 


2,2  2  ,2 

kQr  (AG) 


(4.24) 


lower  diagonal: 


6  1 


4  ky  <48)2 


(4.25) 


The  solution  vector  u^+}  contains  L*M  elements  where  m  =  1,2,...,M  for  each  i  - 

m,  l 

1,2,..., L.  M  is  the  number  of  acoustic  sensors  in  the  vertical  column,  and  L 
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is  the  number  of  vertical  columns  as  varied  from  port  to  starboard  in  azimuth. 

Solution  vector  represents  the  acoustic  field  at  (j+l)*Ar,  m*Az,  and  £*A0 , 

m ,  x, 

i.e.,  range,  depth,  and  azimuth,  respectively. 


Observing  that  the  two  terms  in  square  brackets  on  the  left-hand  side  of 
equation  (4.5)  are  conjugates  of  the  corresponding  two  terms  in  square  brackets 
on  the  right-hand  side,  equation  (4.5)  may  be  expressed  in  the  matrix  form, 


Aj+1Bj+1uj+1  =  AjBjuj 


(4.26) 


or, 


*  *  i+l  -t 

A  B  uJ  =  A  B  uJ 


(4.27) 


where  A  and  B  are  conjugates  of  A  and  B.  Evaluating  the  A  and  B  matrices 
i  i  f  1 

midway  between  uJ  and  uJ  at  range  r+Ar/2  ensures  that  the  operators  A  and  A 
* 

and  B  and  B  are  conjugates  of  each  other. 
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SECTION  5.  COMPUTER  ALGORITHM 


In  the  previous  section,  it  was  shown  how  Lee,  Saad,  and  Schultz  applied 
finite  difference  techniques  to  develop  a  method  for  numerically  solving  the 
LSS  three-dimensional,  wide  angle  wave  equation.  In  this  section  we  will  dis 
cuss  the  algorithm  that  was  used  to  implement  their  method  into  computer  code 

5.1.  TWO-STEP  PROCEDURE  IN  MATRIX  FORM 


Equation  (4.27)  may  be  solved  in  two  steps,  as  follows.  Writing  equation 
(4.27)  in  the  matrix  form 


*  *  i+1  i 

A  B  uJ  =  ABuJ, 


(5.1) 


i+1  *  i+l  i  i 

and  substituting  vJ  for  B  uJ  ,  and  vJ  for  Bu  ,  we  have 


A*vj+1  =  Avj 


(5.2) 


Now,  solve  as  follows: 

Step  1 

Compute  Av^  and,  using  a  tri-diagonal  solver,  solve  the  system 

A  vj+1  =  Av^  (5.3) 


j+1 

for  the  vector  v 
Step  2 

Now  that  vj+1  is  known,  use  a  tri-diagonal  solver  to  solve  the  system 


*  j+1  j+1 

B  u  =  v 


(5.4) 


for  the  vector  u^+^.  As  previously  stated,  u^+^“  is  the  acoustic  field  at 
all  receivers  in  the  cylindrical  plane  at  the  advanced  range,  r+Ar. 
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5.2  COMPUTATIONS  FOR  STEP  1 

Figures  2  and  3  summarize  the  computations  performed  in  Step  1  for  a  sys 
tern  of  equations  where  L“3  and  M=3. 

5.3  COMPUTATIONS  FOR  STEP  2 

Figure  4  summarizes  the  computations  performed  in  Step  2  for  a  system  of 
equations  where  L=3  and  M-3. 
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Figure  2.  STEP  1;  Left-Hand  Side  of  System  of  Equations 
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Figure  3.  STEP  Is  Right-Hand  Side  of  System  of  Equations 
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SECTION  6.  COMPUTER  MODEL 

The  previous  section  presented  the  two-step  marching  procedure  for  solving 
the  LSS  three-dimensional,  wide  angle  wave  equation  and  computations  required 
in  each  of  the  steps.  In  this  section,  we  will  discuss  the  computer  model, 
F0R3D,  which  implements  the  two-step  procedure.  F0R3D  is  written  in  Fortran 
using  single-precision,  complex  arithmetic,  and  has  been  installed  on  a  VAX- 
11/780  digital  computer.  The  model  was  designed,  in  some  instances  at  the 
expense  of  speed,  such  that  additional  capabilities  may  be  easily  incorporated. 
After  the  model  has  been  sufficiently  developed,  it  is  anticipated  that  special 
versions  of  F0R3D  will  be  generated  to  maximize  computational  speed. 

6.1  CAPABILITIES 

The  model  is  capable  of  predicting  acoustic  propagation  loss  in  range-, 
depth-,  and  azimuthal-dependent  ocean  environments.  However,  in  its  present 
stage  of  development,  the  model  is  only  capable  of  handling  a  flat  bottom. 
Boundary  conditions  are  supplied  by  either  the  user  or  the  model  depending  on 
flags  set  in  the  input  runstream.  Additional  optional  boundary  conditions  will 
be  incorporated  into  the  model  as  development  continues.  Methods  for  treating 
irregular  interfaces,  three-dimensional  sound  speed  profile  arrays,  and  other 
starting  fields  are  also  under  consideration. 

6.2  STRUCTURE  OF  MODEL 

The  hierarchical  structure  of  the  model  is  shown  in  figure  5.  Those  sub¬ 
routines  marked  with  an  asterisk  (*)  are  prepared  by  the  user.  Input  parame¬ 
ters  that  result  in  control  being  transferred  to  user- supplied  subroutines  are 
shown  in  the  figure.  For  example,  if  input  parameter  ISF=1,  then  SFLD3D  trans¬ 
fers  control  to  USFLD3D.  If  ISF=0,  subroutine  SFLD3D  returns  a  Gaussian  start- 
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Figure  5. 


Hierarchical  Structure  of  Computer  Model 


28 


TR  7943 


ing  field  to  the  main  program,  F0R3D.  Note  that  the  hierarchical  order  shown 
in  figure  5  is  not  necessarily  the  order  in  which  the  subroutines  are  executed. 

F0R3D  can  be  compiled,  linked,  and  executed  by  submitting  the  following 
commands : 

$FOR  F0R3D , AMIFD3 , BMIFD3 , HNKL , PRINTP , TWOSTEP , RHS , TRID3D , INDX3D , SVP3D , - 

SFLD3D , SC0N3D , BC0N3D , P0RT3D , STBD3D , USVP3D , USFLD3D , USC0N3D , UBC0N3D , - 
UP0RT3D ,USTBD3D , EXACT 

$LINK  F0R3D , AMIFD3 , BMIFD3 , HNKL , PRINTP , TWOSTEP , RHS , TRID3D , INDX3D , SVP3D , - 

SFLD3D , SC0N3D , BC0N3D , P0RT3D , STBD3D , USVP3D , USFLD3D , USC0N3D , UBC0N3D , - 
UP0RT3D , USTBD3D , EXACT 

$RUN  F0R3D 

6.3  GEOMETRY  OF  PROPAGATION 

To  facilitate  subsequent  discussions  on  the  purpose  of  each  sub-routine 
shown  in  figure  3,  the  geometry  through  which  the  acoustic  wave  propagates  and 
the  data  structures  in  which  that  wave  is  contained  are  presented  in  this 
section . 

Figure  6  illustrates  the  top  view  of  a  cylindrical  wedge  through  which  an 
acoustic  wave  is  propagating  from  the  cylindrical  plane  y  to  the  cylindrical 
plane  x. 


Figure  7  shows  a  three-dimensional  view  of  one  sector  of  the  wedge  illus¬ 
trated  in  figure  6. 

Given  the  acoustic  field  along  the  cylindrical  plane  y  at  the  present 
range  (r),  the  purpose  of  the  model  is  to  predict  the  acoustic  field  at  the 
cylindrical  plane  x  at  the  advanced  range,  r+Ar.  Having  solved  for  the  field 
at  range  r+Ar,  the  model  can  then  solve  for  the  field  at  r+2Ar.  In  this  man¬ 
ner,  the  model  marches  forward,  one  step  at  a  time,  until  maximum  range  has 
been  reached. 
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Figure  6.  Top  View  of  Cylindrical  Wedge 
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6.4  DATA  STRUCTURES  OF  THE  ACOUSTIC  FIELD 

The  outwardly  marching  acoustic  field  along  the  cylindrical  plane  y,  as 
viewed  from  the  source,  is  as  shown  in  figure  8.  Note  that  the  left  subscript 
varies  in  depth  from  top  to  bottom  while  the  right  subscript  varies  in  azimuth 
from  port  to  starboard. 

The  data  structures  in  which  the  acoustic  field  is  stored  are  shown  in 
figure  9. 

Although  the  subscript  N  is  used  in  both  figures  to  point  to  the  bottom, 
F0R3D  will,  depending  on  the  type  of  bottom,  use  the  subscript  N+l  to  point  to 
the  bottom  condition. 

6.5  MAIN  PROGRAM  F0R3D 

The  main  program,  F0R3D,  controls  the  execution  of  the  various  subroutines 
that  make  up  the  model.  Initially,  F0R3D  reads  input  parameters  and  environ¬ 
mental  data  from  input  file,  F0R3D.IN,  and  performs  required  initialization. 
Depending  on  input  parameter  ISF,  F0R3D  then  calls  on  either  subroutine  SFLD3D 
or  subroutine  USFLD3D  to  generate  the  starting  field  that  is  to  be  marched  out 
in  range.  Selected  problem  parameters  are  then  printed  and,  if  requested, 
written  in  an  output  file  for  subsequent  use.  F0R3D  then  calls  on  subroutine 
INDX3D  to  compute  an  index  of  refraction  table  that  is  to  be  used  in  the  compu¬ 
tation  of  the  main  diagonals  of  the  tri-diagonal  matrices  A  and  B. 

After  these  preliminary  procedures  have  been  accomplished,  F0R3D  enters 
the  program's  main  loop  and  continues  to  cycle  in  the  loop  until  the  solution 
has  been  marched  out  to  maximum  range,  as  requested. 


31 


TR  7943 


- SURFACE,  SURY - ► 


PORT 

SIDEWALL, 

PORTY 


Ui.o 

1 

f 

c 

z 

© 

STARBOARD 
SIDEWALL,  STBDY 


ACOUSTIC 
FIELD,  U 


- BOTTOM,  BOTY - 
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At  the  beginning  of  each  pass  through  the  loop,  F0R3D  updates  environmen¬ 
tal  parameters,  as  required,  F0R3D  then  calls  on  subroutine  RHS  to  compute  the 
right-hand  side  of  the  system  of  equations  to  be  solved.  Subroutines  INDX3D, 
AMIFD3,  and  BMIFD3  are  then  called  to  compute  the  index  of  refraction  table  and 
the  diagonals  in  the  A  and  B  matrices,  respectively.  F0R3D  then  calls  on  sub¬ 
routine  TWOSTEP  to  march  the  solution  one  range  step  (Ar)  forward.  The  solu¬ 
tion  returned  by  TWOSTEP  is  then  printed  and  written  in  output  file,  F0R3D.0UT, 
as  requested  by  the  user.  When  the  solution  has  been  marched  out  to  maximum 
range,  the  program  is  terminated.  If  the  range  of  the  solution  has  not  reached 
maximum  range,  F0R3D  returns  control  to  the  top  of  the  main  loop  and  repeats 
the  above  procedures.  The  flow-chart  in  figure  10  summarizes  these  procedures. 

6.6  SUBROUTINES 

A  brief  description  of  each  of  the  subroutines  that  make  up  F0R3D  is  pre¬ 
sented  in  this  section. 

6.6.1  Subroutine  SFLD3D 

If  input  parameter  ISF=0,  main  program  F0R3D  calls  on  subroutine  SFLD3D  to 
generate  a  Gaussian  starting  field  at  range  zero.  The  Gaussian  starting 


,  ,1,13  . 
field  is 


U(M+I )  =  CMPLX(PR, 0 . 0) , 


where 


M  =  (J-1)*N 
J  =  1,2, .. . , NSOL 
I  =  1,2, ... ,N 

N  =  number  of  equispaced  points  in  the  vertical  column 


,  •  •  .  , 


NSOL  =  number  of  columns  in  azimuth 


PR  =  GA 
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ZM  “  depth  of  point  =  I*AZ  (meters) 

AZ  -  depth  increment  (meters) 

ZS  =  source  depth  (meters) 

GW  =  Gaussian  width  =  2/FK 
FK  =  reference  wavenumber  =  2*pi*FRQ/C0 
FRQ  =  frequency  (hz) 

CO  =  reference  sound  speed  (meters/sec) 

GA  =  Gaussian  amplitude  “  (  1/GW)*( 2/FK)**0 . 5 

6.6.2  Subroutine  UFLD3D 

If  input  parameter  ISF-1,  main  program  F0R3D  calls  on  user-written  subrou¬ 
tine  UFLD3D  to  generate  the  starting  field.  UFLD3D  must  load  N*NSOL  values  of 
of  the  complex  starting  field  in  array  U. 

If  ISF=0,  UFLD3D  is  not  called  and  may  be  a  dummy  subroutine. 

If  the  user- supplied  starting  field  is  an  elliptic  solution,  then  the 
starting  field  must  be  divided  by  the  Hankel  function.  The  solution  field 
obtained  must  then  be  multiplied  by  the  Hankel  function  before  computing  trans¬ 
mission  loss.  This  procedure  can  be  accomplished  automatically  by  setting 
input  parameter  IHNK  =  1. 

6.6.3  Subroutine  SVP3D 

When  the  range  (RA)  of  the  next  solution  to  be  obtained  is  equal  to  the 
range  of  the  next  set  of  sound  speed  profiles  (RSVP),  the  next  set  of  profiles 
is  input.  If  input  parameter  KSVP=0,  then  subroutine  SVP3D  is  called  to  read 
the  profiles  from  the  input  runstream.  SVP3D  then  loads  the  six  arrays  shown 
in  figure  11. 

Array  ZLYR  contains  the  maximum  depth  of  each  layer  as  measured  from  the 
surface.  Arrays  RHO  and  BETA  contain  density  and  attenuation,  respectively,  in 
each  layer.  Variable  NLYR  is  the  number  of  layers  in  the  profile.  The  profile 
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in  each  layer  is  stored  in  arrays  ZSVP  and  CSVP.  NSVP  is  the  total  number  of 
points  in  the  profile.  If  an  error  is  detected  while  inputting  the  profiles, 
NSVP  is  reset  to  zero.  Array  IXSVP  contains  indices  that  point  to  the  last 
sound  speed  in  each  layer. 

At  the  present  stage  of  development,  interpolation  is  performed  only  in 
depth.  Profile  changes  in  range  are  abrupt,  and  NSOL  profiles,  one  for  each 
sector  boundary,  must  be  provided  since  interpolation  in  azimuth  is  not  per¬ 
formed,  either.  Linear  interpolation  in  depth  is  performed  in  subroutine 
INDX3D . 

Plans  for  the  future  include  the  development  of  a  sound  speed  profile  pre¬ 
processor  which  will  perform  interpolation  in  all  three  dimensions.  Given  a 
three-dimensional  array  of  sound  speed  profiles,  this  pre-processor  will  gene¬ 
rate  profiles  as  required  by  F0R3D  as  it  marches  through  a  three-dimensional 
wedge. 

6.6.4  Subroutine  USVP3D 

If  input  parameter  KSVP  ^  0,  then  subroutine  USVP3D  will  be  called  to  sup¬ 
ply  an  updated  set  of  sound  speed  profiles.  If  KSVP  is  not  reset  to  zero  by 
USVP3D,  then  USVP3D  will  be  called  at  each  step  in  range.  User-prepared  sub¬ 
routine  USVP3D  must  load  the  six  arrays  with  NSOL  profiles  as  described  in  the 
previous  section.  Variables  NLYR  and  NSVP  must  also  be  loaded  by  the  user. 

Variable  KSVP  may  be  used  in  a  computed  GOTO  statement  to  transfer  control 
within  user  program  USVP3D.  When  the  user  no  longer  needs  USVP3D,  KSVP  must  be 
reset  to  zero  within  USVP3D.  The  last  profile  entered  will  be  used  until  the 
solution  range  is  equal  to  the  next  RSVP.  If  KSVP  is  not  reset  to  zero,  then 
USVP3D  will  be  called  until  the  solution  range  is  equal  to  RSVP,  the  range  of 
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1  =  1,2,.  NLYR 

L  =  1,2 _ NSOL 

J  =  1,2 . NSVP 


NLYR  =  NUMBER  OF  LAYERS 

NSOL  =  NUMBER  OF  SECTOR  BOUNDARIES  (SOLUTIONS) 
NSVP  =  NUMBER  OF  POINTS  IN  SOUND  SPEED  PROFILE 


ZLYR  (l,L) 


NLYR 

* 


MAXIMUM  DEPTH 
OF  LAYER  1 
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_  _  J 
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★ 
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IN 
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★ 


*  EXTENDED  IF  USER  REQUESTS  AN  ARTIFICIAL  ABSORBING  LAYER 


Figure  11.  Sound  Speed  Profile  Arrays 
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the  next  set  of  NSOL  profiles.  With  this  option,  the  user  can  generate  a  new 
set  of  profiles  at  each  range  step.  Sound  speed  profile  values  interpolated  in 
all  three  dimensions  may  be  entered  by  the  user  in  this  manner.  When  the  next 
solution  range  is  equal  to  the  next  RSVP,  either  SVP3D  or  USVP3D,  depending  on 
the  next  KSVP,  is  called  to  input  the  next  set  of  profiles. 

6.6.5  Subroutine  PRINTP 

Subroutine  PRINTP  prints  selected  input  parameters. 

6.6.6  Subroutine  INDX3D 

Subroutine  INDX3D  determines  the  values  of  sound  speed  (CL)  and  attenua¬ 
tion  (B)  to  be  used  at  the  receiver  depth  under  consideration.  The  index  of 

2 

refraction  squared  (n  )  is  then  calculated  where  n  =  C^/CL,  and  anc*  CL  are, 
respectively,  a  reference  sound  speed  and  the  speed  of  sound  at  the  depth  under 
consideration.  Linear  interpolation  in  depth  for  CL  is  performed  as  required. 
The  index  of  refraction  squared  is  then  modified  as  follows, 

n2  =  (C./C.)2  +  i(C./C.)2  *  S/27. 287527, 

0  i  0  i 

where 

B  =  attenuation  in  dB/wavelength ,  and 
27.287527  =  pi*20*ln(e). 

Adding  the  imaginary  term  to  the  index  of  refraction  squared  accounts  for 
11^  2 

attenuation.  5  The  new  value  of  n  is  used  in  subsequent  calculations  of  the 

★ 

main  diagonals  in  the  A  and  A  matrices. 
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6.6.7  Subroutine  RHS 

The  purpose  of  subroutine  RHS  is  to  compute  the  right-hand  side  of  the 
system  of  equations  solved  in  STEP  1  of  the  two-step  procedure.  After  calling 
on  SC0N3D,  BC0N3D,  P0RT3D,  and  STBD3D  for  surface,  bottom,  port,  and  starboard 
boundary  conditions,  RHS  performs  the  computations  summarized  in  figure  3.  The 
result  is  stored  in  the  array  (D).  If  the  boundary  conditions  at  the  advanced 
range  are  known,  then  these  conditions  are  included  in  the  calculation  of  the 
right-hand  side.  Plans  for  the  future  include  handling  cases  where  the  boun¬ 
dary  conditions  at  the  advanced  range  are  unknown  but  may  be  expressed  as  a 
function  of  the  solution  to  be  obtained.  Such  a  case  will  require  modifica¬ 
tions  to  some  of  the  matrix  coefficients  on  the  left-hand  side  of  the  system  of 
equations . 

6.6.8  Subroutine  SC0N3D 

SC0N3D  is  called  upon  by  subroutine  RHS  to  supply  surface  conditions.  If 
input  parameter,  ITYPES=0,  a  pressure  release  surface  condition  is  returned; 
i.e.,  the  arrays  SURY  and  SURX  are  set  to  zero.  If  ITYPES=1,  SC0N3D  calls  on 
USC0N3D . 

6.6.9  Subroutine  USC0N3D 

If  ITYPES=1,  the  user  must  write  subroutine  USC0N3D  to  supply  the  surface 
condition.  The  values  of  the  surface  condition  must  be  stored  in  arrays  SURY, 
shown  in  figure  9,  and  in  SURX,  not  shown.  SURY  is  used  to  store  the  surface 
condition  at  the  present  range,  while  SURX  is  used  to  store  the  surface  condi¬ 
tion  at  the  advanced  range.  Surface  conditions,  u^  ^  at  the  present  range  and 
u^+J  at  the  advanced  range,  as  shown  in  figures  2  and  3,  are  used  in  the  compu- 
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tation  of  and  v^+J,  respectively.  The  subscript  &  is  varied  from  0 

U  j  JC  U  j  X 

through  L+l.  If  ITYPES  f  1,  then  USC0N3D  is  not  called  and  may  be  a  dummy 
subroutine . 

6.6.10  Subroutine  BC0N3D 

BC0N3D  is  called  by  subroutine  RHS  to  supply  bottom  conditions.  If  input 
parameter,  ITYPEB  =  0,  the  bottom  is  set  to  zero,  i.e.,  arrays  BOTY  and  BOTX 
are  set  to  zero.  If  ITYPEB  =  1,  then  BC0N3D  calls  on  user  written  subroutine 
UBC0N3D  to  supply  the  bottom  condition.  If  ITYPEB  =  3,  then  BC0N3D  supplies  an 
artificially  absorbing  flat  bottom. 

6.6.11  Subroutine  UBC0N3D 

User-written  subroutine  UBC0N3D  is  called  upon  to  supply  the  bottom  condi¬ 
tion  whenever  input  parameter  ITYPEB  is  set  to  1.  UBC0N3D  must  load  arrays 
BOTY,  shown  in  figures  8  and  9,  and  BOTX  with  bottom  conditions  u^  £  and  u^+J, 
at  the  present  and  advanced  ranges,  respectively.  Bottom  conditions  u^  ^  and 
u^+J  are  used  in  the  computation  of  v^  and  v^  ^ ,  respectively,  as  shown  in 
figures  2  and  3.  Note  that,  in  figures  8  and  9,  the  subscript  l  is  varied  from 
0  through  L+l.  If  ITYPEB  i  1,  then  UBCON3D  is  not  executed  and  may  be  a  dummy 
subroutine . 

6.6.12  Subroutine  P0RT3D 

Subroutine  RHS  calls  on  subroutine  PORT3D  to  supply  the  boundary  condi¬ 
tions  at  the  port  sidewall.  If  input  parameter  ITYPPW  =  0,  PORT3D  sets  the 
arrays  PORTY,  shown  in  figures  8  and  9,  and  PORTX  equal  to  zero.  If  ITYPPW  = 

1,  user-written  subroutine  UP0RT3D  is  called  by  PORT3D.  Future  plans  include 
additional  boundary  condition  options. 
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6.6.13  Subroutine  UP0RT3D 

User-supplied  subroutine  UP0RT3D  is  called  by  P0RT3D  to  supply  port  side- 
wall  boundary  conditions  whenever  input  parameter  ITYPPW  =  1.  UP0RT3D  must 
load  arrays  PORTY  and  PORTX  with  the  wall  conditions  u“l  Q  and  at  the  pres¬ 

ent  and  advanced  ranges,  respectively.  Wall  conditions  u^  ^  and  u^  J  are  used 
in  the  computation  of  v~|  ^  and  v“|+^  as  shown  in  figures  2  and  3.  Note  in  fig¬ 
ures  8  and  9  that  the  subscript  i  is  varied  from  1  through  N-l.  If  ITYPPW  f  1, 
then  UP0RT3D  is  not  called  and  may  be  a  dummy  subroutine. 


6.6.14  Subroutine  STBD3D 


Subroutine  RHS  calls  on  subroutine  STBD3D  to  supply  the  boundary  condi¬ 
tions  at  the  starboard  sidewall.  If  input  parameter  ITYPSW  =  0,  STBD3D  sets 
the  arrays  STBDY,  shown  in  figures  8  and  9,  and  STBDX  equal  to  zero.  If  ITYPSW 
=  1,  user-written  subroutine  USTBD3D  is  called  by  STBD3D.  Future  plans  include 
additional  boundary  condition  options. 

6.6.15  Subroutine  USTBD3D 


User- supplied  subroutine  USTBD3D  is  called  by  STBD3D  to  supply  starboard 

sidewall  boundary  conditions  whenever  input  parameter  ITYPSW  =  1.  USTBD3D  must 

load  arrays  STBDY  and  STBDX  with  the  wall  conditions  u^  .  -  and  -  at  the 

present  and  advanced  ranges,  respectively.  Wall  conditions  u^  and  u^ 

are  used  in  the  computation  of  v?  T  ,,  and  -  as  shown  in  figures  2  and  3. 

r  i,L+l  i,L+l 

Note  in  figures  8  and  9  that  the  subscript  i  is  varied  from  1  through  N-l.  If 
ITYPSW  f  1,  then  USTBD3D  is  not  called  and  may  be  a  dummy  subroutine. 
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6.6.16  Subroutine  AMIFD3 

Subroutine  AMIFD3  computes  the  upper,  lower,  and  main  diagonals  of  the  A 
matrix.  Computations  performed  are  summarized  in  figure  2. 

6.6.17  Subroutine  BMIFD3 

Subroutine  BMIFD3  computes  the  upper,  lower,  and  main  diagonals  of  the  B 
matrix.  Computations  performed  are  summarized  in  figure  4. 

6.6.18  Subroutine  TWOSTEP 

Subroutine  TWOSTEP  marches  the  acoustic  field  one  step  forward  by  calling 
on  the  tri-diagonal  solver,  TRID3D ,  twice.  Given  the  A  matrix  and  the  right- 
hand  side  of  the  system  of  equations  summarized  in  figures  2  and  3,  the  purpose 
of  the  first  call  to  TRID3D  is  to  compute  the  vector  v^+1.  Once  the  vector 
v^+^  has  been  obtained,  TWOSTEP  reorders  the  vector  as  shown  on  the  right-hand 
side  of  the  system  of  equations  shown  in  figure  4.  Given  the  B  matrix  and  the 
reordered  vector  v^+\  the  purpose  of  the  second  call  to  TRID3D  is  to  compute 
the  solution  field  u^+^  at  the  next  step  in  range.  TWOSTEP  reorders  the  solu¬ 
tion  field  uj+1  returned  by  TRID3D  as  shown  in  figure  9. 


6.6.19  Subroutine  TRID3D 

Subroutine  TRID3D,  a  modified  version  of  subroutine  TRIDAG  presented  in 
reference  15,  solves  a  system  of  N  linear  simultaneous  equations  having  a  tri¬ 
diagonal  coefficient  matrix. 

6.6.20  Complex  Function  HNKL 

HNKL  computes  the  Hankel  function.  The  algorithm  for  computing  the  Hankel 
function  is  given  in  reference  16. 
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6.6.21  Subrout ine  EXACT 

Subroutine  EXACT  is  used  for  test  purposes  only.  If  the  exact  solution  is 
known,  subroutine  EXACT  may  be  coded  to  compute  and  print  items  such  as  rela¬ 
tive  error  between  the  computed  and  exact  solutions.  If  the  exact  solution  is 
unknown,  subroutine  EXACT  should  be  a  dummy  subroutine. 

6.6.22  Common  Data  File 

In  order  to  facilitate  further  development  of  the  model,  most  variables 
and  arrays  have  been  declared  in  the  common  file,  COMMON. FOR.  COMMON. FOR  is 
included  in  each  subroutine  as  required. 
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SECTION  7.  USER'S  GUIDE 


The  next  two  sections  describe  input  and  output  formats  in  detail.  A  test 
example  showing  a  sample  input  runstream  and  user-written  subroutines  can  be 
found  later  in  this  report. 

7.1  INPUT  FORMAT 


Prior  to  executing  the  F0R3D  model,  the  input  runstream  containing  problem 
parameters  must  be  stored  in  file  F0R3D.IN.  File  F0R3D.IN  is  assigned  to  FOR¬ 
TRAN  unit  number  NIU  in  the  main  program.  If  the  user  prefers  to  input  problem 
parameters  on  cards,  then  parameter  NIU  should  be  equated  to  the  card  reader 
unit  number,  and  the  statement  which  assigns  file  F0R3D.IN  should  be  removed 
from  the  main  program.  In  either  case,  the  input  runstream  is  prepared  in  free 
format  as  follows. 


7.1.1  Input  Runstream 


Card  Parameters 


1 

2 

3 

4 

5 

6 


N 

N+l 


TITLE 

FRQ , ZS , CO , ISF , RA , ZA , N , IHNK , ITYPES , ITYPEB , ITYPPW , ITYPSW , FLDW , NSEC 
RMAX , DR , WDR , WDZ , PDR , PDZ , I SFLD , ISVP , IBOT 
Ul,U2,U3,U4,U5,U6 
Rl , Zl  ** 

R2  j  Z2  *  BOTTOM  PROFILE 

**  RANGE,  WATER  DEPTH  (METERS) 

★ 

** 

-1,-1 


N+2 

RSVP 

* 

N+3 

KSVP 

* 

************************************** 

* 

N+4 

NLYR 

* 

* 

**************************** 

* 

** 

N+5 

ZLYR( I , L)  , RHO( I ,L) , BETA( I , L) 

* 

** 

*** 

NOTE 

N+6 

zsvp( 1 ,l) ,CSVP( 1 ,l) 

** 

*** 

NOTE 

** 

1 

N+7 

ZSVP(2,L) ,CSVP(2,L) 

**  NOTE 

** 

2 

* 

,  # 

*  3 

* 

* 

N+M 

ZSVP( J,L) ,CSVP(J,L)  ************ 

* 

* 

********************** 

* 

******************************** 
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7.1. 

Card 

1 

2 


Input  parameters  and  notes  are  defined  in  the  next  subsection. 


Input  Parameters 


Units  are  in  meters  and  meters/sec  except  as  noted  below. 

Parameter 

TITLE  =  Comment  card.  80  characters  maximum. 

FRQ  =  Frequency  (HZ) 

ZS  =  Source  depth 

CO  =  Reference  sound  speed.  If  CO  =  0.0,  CO  is  set  to  average 
speed  in  first  layer. 

ISF  =  Starting  field  flag.  0  -  Gaussian.  1  -  User  field.  If  ISF  = 
0,  RA  is  set  to  zero. 

RA  -  Horizontal  range  from  source  to  starting  field.  RA  is  set  to 
0.0  if  ISF  =0. 

ZA  *  Depth  of  starting  field  at  range  RA.  If  ZA  =  0.0,  ZA  is  set 
to  max  depth  of  bottom  layer  in  first  profile.  If  ITYPEB  =  2 
or  3  and  ZA  =  0.0,  ZA  is  set  to  (4/3)*max  depth  of  bottom 
layer.  If  ITYPEB  =  2  or  3  and  ZA  not  zero,  the  artificial 
bottom  layer  is  extended  to  ZA  meters  provided  that  ZA  is 
greater  than  or  equal  to  max.  depth  of  bottom  layer  in  first 
profile . 

N  =  Number  of  equi-spaced  receivers  in  starting  field.  If  N  =  0, 

N  is  set  so  that  the  receiver  depth  increment  is  less  than  or 
equal  to  1/10  wavelength.  If  N  is  greater  than  MXN,  N  is  set 
to  MXN. 

IHNK  =  Hankel  function  flag.  IHNK  =  0,  don't  use  Hankel  function. 

IHNK  =  1,  divide  starting  field  by  Hankel  function,  then  mul¬ 
tiply  the  solution  field  by  Hankel  function  before  computing 
propagation  loss.  If  starting  field  is  Gaussian,  IHNK  should 
be  set  to  0.  If  starting  field  is  elliptic,  IHNK  should  be 
set  to  1. 

ITYPES  -  Type  of  surface 

0  -  Pressure  release.  SC0N3D  sets  SURY  and  SURX  =  0.0 

1  -  User  supplies  surface  condition.  See  subroutine  USC0N3D. 

2  -  Not  implemented  yet. 

ITYPEB  =  Type  of  bottom 

0  -  Pressure  release.  BC0N3D  sets  BOTY  and  BOTX  =  0.0 

1  -  User  supplies  bottom  condition.  See  subroutine  UBC0N3D. 

2  -  Not  implemented  yet. 

3  -  Absorbing  layer  introduced  -  flat  bottom 

4  -  Rigid  bottom  condition  -  Not  implemented  yet 
ITYPPW  =  Type  of  port  sidewall  boundary  condition 

0  -  Field  along  port  sidewall  is  set  to  0.0 

1  -  User  supplied.  See  subroutine  UPORT3D. 

2  -  Not  implemented  yet. 

ITYPSW  =  Type  of  starboard  sidewall  boundary  condition 

0  -  Field  along  starboard  sidewall  is  set  to  0.0 

1  -  User  supplied.  See  subroutine  USTBD3D. 

2  -  Not  implemented  yet. 
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Card 

3 


4 

5 

N+l 

N+2 

N+3 

N+4 

N+5 


Parameter 

FLDW  =  Width  of  field  in  degrees. 

NSE  =  Number  of  sectors  in  field. 

RMAX  =  Maximum  range  of  solution 

DR  =  Range  step.  If  DR  =  0,  DR  is  set  to  1/10  wavelength.  If 

bottom  of  problem  is  not  flat,  DR  is  recomputed  so  that  max 
depth  is  either  incremented  or  decremented  by  DZ.  Solution  is 
computed  every  DR  meters. 

WDR  =  Range  step  at  which  solution  is  written  on  disk.  If  WDR  not 

0,  an  output  disk  file  is  assigned.  WDR  is  rounded  to  nearest 
DR. 

WDZ  =  Depth  increment  at  which  solution  is  written  on  disk.  Rounded 
to  nearest  DZ. 

PDR  =  Range  step  at  which  solution  is  printed.  Rounded  to  nearest 
DR. 

PDZ  =  Depth  increment  at  which  solution  is  printed.  Rounded  to 
nearest  DZ. 

ISFLD  =  0  -  don't  print  starting  field. 

1  -  print  starting  field. 

ISVP  =  0  -  don't  print  sound  velocity  profile. 

1  -  print  sound  velocity  profile. 

IBOT  =  0  -  don't  print  bottom  depths. 

1  -  print  bottom  depths. 

U1-U6  =  User  variables  -  real,  single  precision. 

R1,Z1  =  Bottom  profile.  First  range  and  depth  of  water. 

R2,Z2  =  etc. 


-1,-1  =  Marks  the  end  of  the  bottom  profile. 

RSVP  =  Range  of  SVP . 

See  Note  1. 

KSVP  =  SVP  flag. 

=  0  -  SVP  in  input  runstream. 

=  Not  zero.  Profile  (cards  N+4  thru  N+M)  is  supplied  by  user. 
User  writes  subroutine  USVP3D.  KSVP  may  be  used  in  computed 

GOTO  statement  to  transfer  control  in  USVP3D. 

\ 

NLYR  =  Number  of  layers.  If  ITYPEB  *  2  or  3,  program  inserts  an 
artificial  layer  and  increments  NLYR  by  1.  See  Note  2. 

ZLYR(I,L)  =  Max  depth  of  layer  I  in  profile. 

RH0(I,L)  =  Density  in  layer  I  (g/cm**3). 

BETA(I,L)  =  Attenuation  in  layer  I  (dB/wavelength) .  If  BETA(I,L)  is 
negative,  attenuation  is  computed.  See  Note  3. 


47 

* 


TR  7943 


Card  Parameter 

N+6  ZSVP( 1 ,  L)  =  Depth  to  top  of  layer  I 

CSVP(1,L)  =  Speed  of  sound  at  top  of  layer  I 


N+M  ZSVP ( J , L)  =  Depth  to  bottom  of  layer  I 

CSVP(J,L)  ■  Speed  of  sound  at  bottom  of  layer  I.  If  only  one  SVP  input¬ 
ted,  it  is  used  through  entire  problem.  If  more  than  one  SVP 
inputted,  last  SVP  is  used  through  remainder  of  problem. 

Note  1  Repeat  cards  N+2  through  N+M  for  each  set  of  profiles  at  range  RSVP. 

A  set  of  profiles  consists  of  NSOL  (NSEC+1)  profiles  located  on  sector 
boundaries  and  equi-spaced  from  port  to  starboard  in  azimuth.  NSEC  is 
the  number  of  sectors  in  the  field  of  interest.  NSEC  adjacent  sectors 
have  NSOL  unique  boundaries. 

Note  2  Repeat  cards  N+4  through  N+M  for  each  profile  in  the  set.  At  range 
zero  all  profiles  are  the  same.  Beyond  range  zero,  L  profiles  are 
entered  from  port  to  starboard  as  viewed  from  range  RO.  L  is  varied 
from  1  through  NSOL.  If  all  profiles  are  the  same,  enter  -1,-1  after 
the  last  ZSVP  and  CSVP  of  the  first  profile.  This  will  cause  the 
first  profile  to  be  entered  repeatedly  until  NSOL  profiles  have  been 
entered . 

Note  3  Repeat  cards  N+5  through  N+M  for  each  layer  I  in  the  profile.  J  is 
the  number  of  points  in  layer  I. 


7 . 2  OUTPUT  FORMAT 

Output  data  generated  by  the  F0R3D  model  are  written  on  disk  in  file 
F0R3D . OUT.  File  F0R3D . OUT  is  assigned  to  Fortran  unit  number  NOU  in  the  main 
program.  The  data  are  written  with  Fortran  unformatted  WRITE  statements  as 
explained  below. 

The  values  of  the  variables  and  arrays  listed  in  the  following  WRITE 
statement  are  written  once  at  the  beginning  of  the  program  as  follows: 

WRITE (NOU)FRQ , ZS , CO , I SF , RA , ZA , N , IHNK , ITYPES , ITYPEB , ITYPPW , ITYPSW , FLDW , 
NSEC, NSOL, RMAX, DR, WDR,WDZ,DZ,NLYR,ZLYR,RHO, BETA, U1,U2,U3,U4,U5,U6 

These  variables  and  arrays  are  defined  in  the  previous  section. 
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The  values  of  the  variables  and  arrays  listed  in  the  following  WRITE 
statement  are  written  at  RO  and  at  each  WDR  (write  range  increment)  thereafter. 


WRITE ( NOU ) ANG, NZ ,RA, WDZ , (U(M+l) , I=IWZ ,N, IWZ) 


where  M 
J 
N 

ANG 

NZ 

RA 

WDZ 

WDZ 

DZ 

U 


(J-1)*N 


Solution  index.  J  varies  from  1  to  NSOL.  Port  to  starboard. 
Number  of  receivers  in  vertical  column. 

Relative  Bearing  of  solution  in  degrees.  Port  is  negative. 
Number  of  receivers  written  on  disk. 

Range  in  meters. 

Depth  increment  of  receivers  written  on  disk  in  meters. 


IWZ*DZ 


Depth  increment  of  receivers  in  vertical  column. 
Solution  array.  Complex. 


Printed  output  is  directed  to  printer  unit  number  NPU. 


7 . 3  PLOT  PROGRAM 


Since  plotting  facilities  vary  from 
ed  that  plotting  capabilities  should  not 
A  separate  plot  program  has  been  written 


one  activity  to  another,  it  was  decid- 
be  incorporated  into  the  F0R3D  model, 
to  plot  the  output  of  the  model. 
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SECTION  8.  TEST  PROBLEM 

In  this  section  a  test  problem  previously  reported  by  Siegmann,  Lee,  and 
Botseas^  is  included  to  demonstrate  the  accuracy  and  capabilities  of  the 
model.  The  input  runstream  and  the  subroutines  required  to  supply  the  starting 
field,  sound  speed  profiles,  bottom,  surface,  and  sidewall  boundary  conditions 
are  also  included. 

8.1  EXACT  SOLUTION 

The  exact  solution  selected  to  demonstrate  the  accuracy  of  the  model  is 

2  i(M.9  -  R.(r)  +  R . (r  ) ) 

u(r,0,z)  =  £  e  ^  ^  ^  Z,(z)5  (8.1) 

j=i  J 


where 


and 


Z  .  (z)  =  a .  sin(y  .  k_  z) , 
J  J  J  0 


U  .  =  (N .  -  1/2)  TT/Ck-H) 
J  J  0 


R.(r) 

J 


r  1  2  i 
1  2  uj  k0  r  ‘ 


M. 


2ko  r 


[1 


l 

4 


2,-1 

yi] 


In  equations  (8.1)  through  (8.4), 

=  an  integer  azimuthal  mode  number, 
a j  =  a  relative  mode  amplitude, 
r  =  range , 
z  =  depth, 

0  =  azimuthal  angle, 


(8.2) 

(8.3) 


(8.4) 
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kQ  =  2*pi*f/C0, 

=  initial  range, 

=  integer  vertical  mode  number, 

H  =  depth  of  bottom, 

Cq  =  reference  sound  speed,  and 
f  =  source  frequency. 

Input  parameters  selected  for  the  test  case  are  as  follows: 

Mx  =  24,  M2  =  36,  Nx  =  3,  N2  =  5,  a 1  -  1,  a 2  *  1, 

Source  frequency  =  50  hz, 

Source  depth  =  10  meters, 

Initial  range  of  starting  field  =  500  meters, 

Depth  of  water  *  100  meters, 

Sound  speed  in  water  =  1500  m/s , 

Reference  sound  speed  «  1500  m/s , 

Density  in  water  =  1  gm/ cm**3 , 

Attenuation  in  water  =  0  db/wavelength , 

Depth  increment  =  .5  meters, 

Range  increment  -  2  meters, 

Maximum  range  =  1000  meters, 

Width  of  field  in  azimuth  *  10  degrees,  and 
Azimuthal  increment  =  0.1  degree. 


Three-dimensional  plots  of  exact  propagation  loss  solutions  at  20  and  60 
meters  in  depth  are  shown  in  figures  12  and  13,  respectively.  Propagation  loss 
is  computed  as  follows: 

PL  =  -20*LOG10( | u | )  +  10*LOG10(r) . 
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Figure  12.  Exact  Solution  at  20  Meters  in  Depth 
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Figure  13.  Exact  Solution  at  60  Meters  in  Depth 
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As  can  be  seen  in  the  figures,  the  exact  solutions  to  this  test  case  are  spi 
rals  that  radiate  from  the  center  and  vary  in  all  three  dimensions. 


8.2  COMPARISON  OF  RESULTS 


Solutions  obtained  with  F0R3D  along  relative  bearing  lines  -5°,  0°,  and 
+5°  at  receiver  depths  20  and  60  meters  and  over  the  range  interval  500-1000 
meters  as  shown  in  figures  12  and  13  are  compared  with  corresponding  exact 
solutions  in  figure  14. 


Solutions  obtained  with  F0R3D  are  compared  with  each  other  in  figure  15. 


8 . 3  INPUT  RUN STREAM 


The  input  runstream  for  the  test  case  is  listed  below. 


Test  Case  3 

50  10  1500  1  500  100  200  0  0  1  1  1  10  100 
1000  2  2  10  5000  25  0  0  0 
3  5  24  36  1  1 
0  100 
10000  100 
-1,-1 
0 
1 


8.4  USER-SUPPLIED  SUBROUTINES 


The  subroutines  which  were  supplied  for  this  test  case  have  been  merged 
into  one  file,  TEST3.F0R  listed  below.  The  listing  of  TEST3.F0R  is  followed  by 
the  command  file  TEST3.COM  that  compiled,  linked,  and  executed  FOR3D. 
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Figure  14.  Comparison  of  F0R3D  and  Exact  Solutions 
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Figure  15.  Comparison  of  F0R3D  Solutions 
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SUBROUTINE  EXACT 

C  *************************************************** 

C  ***  EXACT  SOLUTION  FOR  SIEGMANN  TEST  CASE 

C  *************************************************** 

INCLUDE  ' COMMON . FOR ' 

SMU1 =(U1-.5)*PI/(ZA*XK0) 

TMP= ( 1 . 0- . 25*SMU1 *SMU I ) 

R 1 R0= ( . 5  *SMU1 *SMU1 *XK0  *R0- (U3  *U3 ) / ( 2 . 0*XK0*R0 ) ) /TMP 
R1RA= ( . 5*SMU1 *SMU1 *XK0 *RA- ( U3*U3 ) /( 2 . 0*XK0*RA ; ) /TMP 
SMU2= ( U2~ . 5 ) *PI / ( ZA*XK0 ) 

TMP= ( 1  . 0- . 2  5  *SMU2  *SMU2 ) 

R2R0=i  (  . 5*SMU2*SMU2  *XK0  *R0  -  (  U4  *U4  )  /  (  2 . 0*XK0  *R0  )  )  /TMP 
R2RA- ( . 5  *SMU2  *SMU2  *XK0  *RA- ( U4  *  U4 ) / ( 2 . 0  *XK0  *RA ) ) /TMP 
SZ1 =U5*S I N ( SMU 1 *XK0  *Z I ) 

SZ2  =  U6*SIN( SMU2  *XK0  *Z I ) 

ANGR=ANG*PI / 180.0 

UB=SZ 1 *CEXP ( CMPLX (0.0, U3  * ANGR-R 1 RA  +  R1 R0 ) ) * 

C  SZ2*CEXP ( CMPLX ( 0 . 0 , U4  *ANGR-R2RA*R2R0 ) ) 

PL  =  -20.0*  ALOG 1 0 (CABS(UB) ) +1 0 . 0  *ALOG 1 0 ( RA ) 

UC=UA-UB 

WRITE ( 6 , 120)  PL , UB 

I F ( ANGR . EQ . 0 ) WRITE ( 6 , * ) RA , Z I , PL , UB 
120  FORMAT (20X,3X,F10.3,3X, ' ( ' ,E12.5,2X,E12.5, '  )') 

RSLERR=CABS ( UC/UB ) 

WR ITE ( 6 , 150)  UC , RELERR 

150  FORMAT ( 3  6X ,  ' ( ' , E 1 2 . 5 , 2X , E 1 2 . 5 , '  )',?X,E12.5) 

RETURN 
END 

C 
C 

SUBROUTINE  USVP3D 

Q  ****************************************************************** 

C  ***  SOUND  VELOCITY  PROFILE  SUBROUTINE 

C  ****************************************************************** 

INCLUDE  'COMMON. FOR' 

GO  TO  (100,200,300,400)  , KSVP 
NSVP=0 
RETURN 
C 

100  CONTINUE 

C  ***  SIEGMANN  TEST  CASE 

NSVP- 2 
NLYR= 1 

DO  75  L= 1 , NSOL 
ZLYR ( 1 ,L) =100.0 
RHO ( 1  ,L)  =  1  .0 
BETA ( 1 , L) =0 . 0 
ZSVPd  ,  L )  =  0 . 0 
CSVP( 1 , L ) = 1 500 . 0 
ZSVP ( 2 , L ) = 1 00 . 0 
CSVP ( 2 , L ) - 1 500 . 0 
1XSVP( 1 , L)=NSVP 
75  CONTINUE 

RETURN 
C 

200  CONTINUE 

C  ***  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 
C 

300  CONTINUE 


;  ******* 
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C  ***  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 
C 

400  CONTINUE 

:  ***  USER  INSERTS  CODE  HERE  IF  DESIRED 

RETURN 
END 


C 

C 

C 


100 


c 

c 

c 

c 

c 


c 

c 

100 

c 

c 

c 

c 

200 


SUBROUTINE  USFLD3D 


***  EXACT  STARTING  FIELD 


INCLUDE  'COMMON. FOR' 

SMU 1— (Ul— .5)*PI/( ZA*XK0 ) 

TMP= ( 1 . 0-. 25*SMU1 *SMU 1 ) 

R1 R0= ( . 5*SMU1 *SMU1 *XK0 *R0- ( U3 *U3 ) / ( 2 . 0*XK0 *R0 ) ) /TMP 
R1RA-( . 5*SMU1 *SMU 1 *XK0 *RA- ( U3 *U3 ) / ( 2 . 0 *XK0 *RA  ) /TMP 
SMU2  =  (U2-.5)*PI/( ZA*XK0 ) 

TMP= ( 1 . 0- . 2  5*SMU2  *SMU2 ) 

R2R0=  (  .  5*SMU2  *SMtJ2  *XK0  *R0-  (U4*U4)/(2. 0*XK0  *R0  )  )  /TMP 
R2RA- ( . 5*SMU2*SMU2*XK0*RA- (U4*U4 ) / ( 2 . 0*XK0*RA) ) /TMP 
DO  100  J=1 , NSOL 
M= ( J- 1 ) *N 


ANG--FLDW/2 .0*PI/180.0+( (J-1 ) *PH  I ) 
DO  100  1=1 ,N 
Z I = I *DZ 


SZ1 =U5  *S I N ( SMU1 *XK0  *Z I ) 

SZ2=U6*SIN(SMU2*XK0*ZI ) 

U( I+M) =SZ1 *CEXP(CMPLX( 0 . 0 ,U3*ANG-R1 RA+R1 R0 ) ) + 
C  U5*SZ2  *CEXP ( CMPLX (0.0, U4 *ANG-R2RA^R2R0 ) ) 

CONTINUE 
RETURN 
END 


SUBROUTINE  UBCON3D 

***  EXACT  BOTTOM  CONDITION  SUBROUTINE 

********»*u***H»»*****«*u****A*«a****»******»|r*»**»(*t(|*jjH|t 

INCLUDE  'COMMON. FOR’ 

I F ( THETA )  100,200,300 

***  THETA  LESS  THAN  0.0.  BOTTOM  SLOPES  UP. 

CONTINUE 
BOTY  =  U ( N ) 

BOTX= . 

RETURN 


0.0.  BOTTOM  IS  FLAT, 


**■*  theta 
CONTINUE 
SMU I  =  ( U 1  - . 5 ) *p ; / ( ZA*XK0 ) 

TMP  --  ( 1  .  0  -  .  2  5  *  SMU  !  '  SM!  i  1  ) 

R1  R0=  (  .5  *  SMU  1  *SMU  1  *XK.0*R0  t  J  3*U  3  )/(  2 . 0  *  XIC0  *R0  )  )  /TMP 
R1  KAX  (  .  b*SMU1  *SMU1  *XKJ*KA  ( <)  i  *  \l.i  )  /  (  2 . 0  *XKG  *HA )  ) /TMP 

^Y7';5^MU1‘SMU'  *XK0M  Pa'or>  -  (U3*U3  )  /(  2 . 0  *XK0 *  (  RA-DR)  )  )  /TMP 
SMU  2  = ( U2 - . 5 ; *  p I / ( Z A  *  X K 0 )  ' 

TMP=  (  I  .  0-  .  2 S *SMU2 XSML'2  ) 

R2R0=(  .  5*SMU2*SMU2*XK0*Rij-  (U4*U4  )/(  2  0*X-0'R0)  )  /'"MP 
R2RAX  (  .  5*SMU?*SMU2*XK*J*RA-  (  U4*'j-1  )/(2 .0*XK0*RA)  )  /TMP 
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R2RAY=( .5*SMU2*SMU2*XK0*(RA-DR)-(U4*U4)/(2.0*XK0*(RA-DR) ) )/TMP 
ZI=(N+1 ) *DZ 
DO  250  J=1  , NSOL-t-2 

ANG=-FoDW/2 .0*PI/180.0+( ( J- 1 ) *PHI ) 

SZ1=U5*SIN( SMU 1 *XK0*ZI ) 

SZ2=U6*SIN(SMU2*XK0*ZI ) 

BOTX ( J ) =SZ1 *CEXP{ CMPLX (0.0 , U3*ANG-R1 RAX+R1 R0 ) ) + 

C  U5  *SZ2  *CEXP ( CMPLX (0.0, U4  *ANG-R2RAX  +  R2R0 ) ) 

BOTY ( J ) =SZ1 *CEXP ( CMPLX (0.0 , U3  *ANG-R1 RAY+R1 R0 ) )  + 

C  U5  *SZ2  *CEXP ( CMPLX (0.0, U4  *ANG-R2RAY  +  R2R0 ) ) 

250  CONTINUE 
RETURN 

***  THETA  GREATER  THAN  0.0,  BOTTOM  SLOPES  DOWN. 

CONTINUE 

BOTY= . 

BOTX= . 

RETURN 
END 


SUBROUTINE  UPORT3D 

***  EXACT  PORT  BOUNDARY  CONDITION 

****************************************  ********if**irir************* 
INCLUDE  'COMMON. FOR' 

SMU1  =  ( U 1  - . 5 ) *P I / ( ZA*XK0 ) 

TMP=(1 . 0- . 25*SMU1 *SMU 1 ) 

R1 R0= ( . 5*SMU1  *SMU  1  *XK0  *R0 - ( U 3*U3 ) / ( 2 . 0*XK0  *R0 ) ) /TMP 
R 1  RAX  = ( . 5*SMU1 *SMU1  *XK0 * RA- ( U 3 *U3 ) / ( 2 . 0*XK0*RA) ) /TMP 
R1 RAY= ( . 5*SMU1 *SMU1 *XK0 * ( RA-DR ) - ( U3 *U3 ) / ( 2 . 0 *XK0 * ( RA-DR ) ) ) /TMP 
SMU2= (U2-.5)*PI/( ZA*XK0 ) 

TMP=(1  . 0- . 2  5*SMU2  *SMU2 ) 

R2R0=( . 5*SMU2*SMU2*XK0*R0- ( U4*U4 ) /( 2 . 0*XK0*R0 ) ) /TMP 
R2RAX= ( .  5*SMU2*5MU2*XK0*RA- (U4*U4 )/ ( 2 . 0*XK0*RA) ) /TMP 
R2RAY=(  . 5* SMU 2  *SMU2  *XK0  * ( RA-DR )-(U4*U4)/(2.0  *XK0  * ( RA-DR )  ) ) /TMP 
ANG= -FLDW/2 . 0  *PI/ 1 80 . 0-PHI 
DO  250  1=1 ,N 
Z I = I *DZ 

SZ1 =U5*S IN ( SMU 1 *XK0*ZI ) 

SZ2=U6*SIN( SMU2  *XK0  *Z I ) 

PORTX  (  I  )  -SZ 1  *CF.XP  (  CMPLX  (  C  .0  ,  U3*ANG~R1  RAX  +  R1  R0  )  )  + 

C  U5  *SZ2 *CEXP( CMPLX ( 0 . 0 , U4  *ANG -R2RAX  +  R2R0 ) ) 

PORTY ( I ) =SZ1 *CEXP( CMPLX ( 0 . 0 , U3*ANG -R 1  RAY *R 1 R0 ) )  + 

C  U5*SZ2*CEXP(CMPLX(0.0,U4*AtJG-R2RAY+R2R0)  ) 

250  CONTINUE 

RETURN 
END 
C 
C 

SUBROUTINE  USTBD3D 

C  *************  A  *********************************************  ******* 

C  ***  EXACT  STARBOARD  BOUNDARY  CONDITION 

C  ****************************  *  ******************************  ******* 

INCLUDE  'COMMON. FOR' 

SMU 1  =  ( U 1  - . 5 ) *P I / ( ZA*XKG ) 

TMP= ( 1 . 0  - . 2  5  *  SMU 1  *  SMU 1 ) 

R1 R0= ( . 5* SMU I  *  SMU 1 *XK  0  *  RO - ( U  3  *  U  3 )/( 2 . 0*XK0*R0 ) ) /TMP 

R1 RAX- ( . 5* SMU 1  * SMU i *XKU  * RA- ( U3 * U  3 ) / ( 2 . 0 *XK0  *RA ) ) /TMP 

R 1  RAY  -  (  .  5*SMU1  *3M'.T  *XK0*  (  RA-DR  J  -  (U3*U  3  )  / «  2 . 0  *XK0*  (  RA-DR )  )  ) /TMP 
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250 


SMU  2  =  ( U  2  - . 5) *PI/ (ZA*XK0 ) 

TMP=(1  .0-.25*SMU2*SMU2)  u4)/(2^QirXK0*R0)  )/TMP 

il^-l  :i^52*SS2*5KO*<Si-SR)-ti4iSi)/c5?5^0M^-DR> )  )/THP 

ANG=FLDW/2.0*PI/180.0+PHI 

DO  250  1=1 ,N 
Z I - I *DZ 

SZ1 =U5*SIN(SMU1 *XK0*ZI ) 


CONTINUE 

RETURN 

END 


61 


TR  7943 


j  COMMAND  FILE  TO  COMPILE,  LINK  AND  EXECUTE  TEST  3 

$SET  VERIFY 
$ ASSIGN  NL  SYS5PRINT 

IZ  JSsSSJ ;S:SBmSF03  .HNKB .  PRINTP .  - 
^SFlKSI  S&3D , PORT3D ,  STBD3D ,  - 
$l“NFOT3d!mJ  FD3  ,  BMIFD3  ,  HKKL ,  PRINTP ,  - 

SVP^D^SFLDSdTsCON^D^ BCON3D~PORT3D, STBD3D, " 

USCON3D, TEST3 

$COPY  TEST 3 . IN  FOR3D.IN 

$RUN  FOR3D 


C 
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SECTION  9.  CONCLUSION 

A  computer  model,  FOR3D,  was  developed  for  implementing  the  Lee-Saad- 
Schultz  method  for  solving  the  LSS  three-dimensional,  wide  angle  wave  equation. 
The  model  is  designed  to  predict  acoustic  propagation  loss  in  range-,  depth-, 
and  azimuthal-dependent  ocean  environments.  The  computational  speed  of  the 
model  is  favorable  since  the  Lee-Saad-Schultz  method  requires  only  solving  two 
tri-diagonal  systems  of  equations  for  each  step  marched  forward  in  range. 

Theory  shows  that  the  Lee-Saad-Schultz  method  is  unconditionally  stable. 

Although  the  model  is  still  under  development,  it  is  presently  capable  of 
accepting  arbitrary  initial,  surface,  bottom,  and  sidewall  boundary  conditions. 
Another  important  feature  of  the  model  is  that  it  has  wide  angle  capability  in 
the  vertical  plane.  Angles  of  propagation  as  wide  as  31  degrees  have  been 
accurately  handled  by  the  model. 

The  model  is  designed,  in  some  instances  at  the  sacrifice  of  speed,  such 
that  future  modifications  and  capabilities  may  be  incorporated  easily.  After 
the  model  has  been  sufficiently  developed,  it  is  anticipated  that  special  ver¬ 
sions  of  the  model  will  be  generated  that  maximize  computational  speed. 

Plans  for  further  development  of  the  model  include  methods  for  treating 
irregular  interfaces,  a  variety  of  optional  starting  fields,  a  variety  of 
optional  boundary  conditions,  and  a  rigorous  treatment  of  a  three-dimensional 
array  of  sound  speed  profiles.  User  contributions  that  will  enhance  the  model 
are  invited. 
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